
// Polynomial template class

// $Id: Poly.hpp 2 2012-05-22 17:55:16Z gerry@gaboury.biz $

#ifndef _POLY_HPP
#define _POLY_HPP

#include <stdexcept>
#include <cmath>
#include <complex>
#include <iostream>
#include <limits>

#include "Vector.hpp"

using namespace std;

template <class T, class X> class Poly
{
protected:
  int n;
  T *C;
public:
  // constructors
  Poly () : n(0), C(0) { }
  Poly ( int order ) : n(order), C(0) {
    if (order < 0)
      throw (range_error("Poly::Poly: order less than zero"));
    else {
      C = new T [n+1];
      for (int k=0; k <= n; ++k)
	C[k] = 0;
    }
  }
  Poly (const T *V, int order) : n(order), C(0) {
    if (order < 0)
      throw (range_error("Poly::Poly: order less than zero"));
    else {
      C = new T [n+1];
      for(int k=0; k <= n; ++k)
	C[k] = V[k];
    }
  }
  Poly ( const vector<T> vec ) : n(vec.size()), C(0) {
    if ( n < 0 )
      throw (range_error("Poly::Poly: order less than zero"));
    else {
      C = new T [n+1];
      for(int k=0; k <= n; ++k)
	C[k] = vec[k];
    }
  }
  // copy constructor
  Poly (const Poly<T,X>& P) : n(P.n) {
    C = new T [n+1];
    for(int k=0; k <= n; ++k) 
      C[k] = P.C[k];
  }
  virtual ~Poly () {
    if (n)
      delete [] C;
  }
  // functions
  int order() const {
    return n;
  }
  void Coeff(T src,int k) {
    if ( k < 0 || k > n )
      throw (range_error("Poly::Coeff - dimension error"));
    C[k] = src;
  }    
  T Coeff(int k) const {
    if (k < 0 || k > n)
      throw (range_error("Poly::Coeff - dimension error"));
    return C[k];
  }
  // derivative
  Poly<T,X> Deriv() {
    if (n < 1) {
      Poly<T,X> P(1);
      return P;
    }
    Poly<T,X> P(n-1);
    for (int k=0; k <= n-1; ++k)
      P.C[k] = C[k+1] * (k+1);
    return P;
  }
  // Indefinite Integral
  Poly<T,X> Integral( T C0=0 ) {
    Poly<T,X> P(n+1);
    for (int k=0; k <= n; ++k)
      P.C[k+1] = C[k] / (k+1);
    P.C[0] = C0;
    return P;
  }
  // shift coefficients (increase the order)
  // This is effectively the same as multiplication by x**m.
  Poly<T,X> Shift(int m) {
    Poly<T,X> P(n+m);
    for (int k=0; k <= n; ++k)
      P.C[k+m] = C[k];
    return P;
  }
  // same as operator ()
  X eval(X x) {
    X tmp = C[n];
    for (int k = n-1; k >=0; --k) {
      tmp *= x;
      tmp += C[k];
    }
    return tmp;
  }
  // operators
  X operator()(X x) {
    X tmp = C[n];
    for (int k = n-1; k >=0; --k) {
      tmp *= x;
      tmp += C[k];
    }
    return tmp;
  }
  T& operator[](int k) {
    if (k < 0 || k > n)
      throw (range_error("Poly::Coeff - dimension error"));
    return C[k];
  }
  Poly<T,X> operator+(Poly<T,X>& P1) {
    int i;
    if (n >= P1.n) {
      Poly<T,X> P2 = *this;
      for (i = 0; i <= P1.n; ++i)
	P2.C[i] += P1.C[i];
      return P2;
    } else {
      Poly<T,X> P2 = P1;
      for (i = 0; i <= n; ++i)
	P2.C[i] += C[i];
      return P2;
    }
  }
  Poly<T,X> operator-(Poly<T,X>& P1) {
    int i;
    if (n >= P1.n) {
      Poly<T,X> P2 = *this;
      for (i = 0; i <= P1.n; ++i)
	P2.C[i] -= P1.C[i];
      return P2;
    } else {
      Poly<T,X> P2 = P1;
      for (i = 0; i <= n; ++i) {
	P2.C[i] *= -1;
	P2.C[i] += C[i];
      }
      return P2;
    }
  }
  Poly<T,X> operator*(Poly<T,X>& P1) {
    Poly<T,X> P2(P1.n + n);
    for (int i=0; i <= n; ++i) {
      for (int j=0; j <= P1.n; ++j)
	P2.C[i+j] += C[i] * P1.C[j];
    }
    return P2;
  }
  // The denominator is a factor of the numerator
  // therefore the division produces a polynomial of
  // known length with no residue.
  // This is used in root finding ... the denominator
  // is either a zero, first or second order polynomial.
  // It is assumed that the highest coefficient of the 
  // denominator is one.
  Poly<T,X> operator/(Poly<T,X>& D) {
    Poly<T,X> P(n - D.n);
    if (P.n <= 0)
      throw (range_error("Poly::Division - numerator order"));
    if (D.n < 0 || D.n > 2) {
      throw (range_error("Poly::Division - denominator order"));
    } else if (!D.n) {
      for (int i=0; i <= n; ++i)
	P.C[i] /= D.C[0];
    } else if (D.n == 1) { // first order denominator
      P.C[n-1] = C[n];
      for (int i=n-2; i >= 0; --i)
	P.C[i] = C[i+1] - D.C[0] * P.C[i+1];
    } else {	// second order denominator
      P.C[n-2] = C[n];
      P.C[n-3] = C[n-1] - D.C[1] * P.C[n-2];
      for (int i=n-4; i >= 0; --i)
	P.C[i] = C[i+2] - D.C[0] * P.C[i+2] - D.C[1] * P.C[i+1];
    }      
    return P;
  }
  /* Not ready for prime time yet.
  //
  // Inverse - reciprocal
  //
  Poly<T,X> Inverse() {
    Vector<T> B();
    // C0
    B.push_back(1.0/C[0]);

    if (P.n <= 0)
      throw (range_error("Poly::Division - numerator order"));
    if (D.n < 0 || D.n > 2) {
      throw (range_error("Poly::Division - denominator order"));
    } else if (!D.n) {
      for (int i=0; i <= n; ++i)
	P.C[i] /= D.C[0];
    } else if (D.n == 1) { // first order denominator
      P.C[n-1] = C[n];
      for (int i=n-2; i >= 0; --i)
	P.C[i] = C[i+1] - D.C[0] * P.C[i+1];
    } else {	// second order denominator
      P.C[n-2] = C[n];
      P.C[n-3] = C[n-1] - D.C[1] * P.C[n-2];
      for (int i=n-4; i >= 0; --i)
	P.C[i] = C[i+2] - D.C[0] * P.C[i+2] - D.C[1] * P.C[i+1];
    }      
    return P;
  }
  */
  void operator+=(const Poly<T,X>& P1) {
    int i;
    if (n >= P1.n) {
      for (i=0; i <= P1.n; ++i)
	C[i] += P1.C[i];
    } else {
      Poly<T,X> P = P1;
      for (i=P1.n; i <= n; ++i)
	P.C[i] += C[i];
      delete [] C;
      n = P.n;
      C = new T [n+1];
      for (i=0; i <= n; ++i)
	C[i] = P.C[i];
    }
  }

  void operator-=(const Poly<T,X>& P1) {
    int i;
    if ( n >= P1.n ) {
      for (i=0; i <= P1.n; ++i)
	C[i] -= P1.C[i];
    } else { // P1.n > n
      Poly<T,X> P ( P1.n );
      for ( i=0; i <= n; ++i )
	P.C[i] = C[i] - P1.C[i];
      for ( i=n+1; i <= P1.n; ++i )
	P.C[i] = -P1.C[i];
      delete [] C;
      n = P.n;
      C = new T [n+1];
      for (i=0; i <= n; ++i)
	C[i] = P.C[i];
    }
  }

  Poly<T,X> operator*(const T c) {
    Poly<T,X> P(n);
    for (int i=0; i <= n; ++i)
      P.C[i] = C[i] * c;
    return P;
  }
  void operator*=(const T c) {
    for (int i=0; i <= n; ++i)
      C[i] *= c;
  }
  virtual Poly<T,X>& operator=(const Poly<T,X>& Ps) {
    if (this == &Ps)
      return *this;	// don't assign to self
    if (n != Ps.n) {
      delete [] C;
      n = Ps.n;
      C = NULL;
    }
    if (C == NULL)
      C = new T [n+1];
    for (int k=0; k <= n; ++k)
      C[k] = Ps.C[k];
    return *this;
  }
  // IO
  void Write(ostream *os) {	// Write
    for (int k=0; k <= n; ++k)
      *os << C[k] << "  ";
    *os << endl;
  }
  inline T NewtonRaphson ( const T& p_init, const T& p_eps, 
			   int n_relax, int n_max )
  {
    T x=p_init, eps1 = p_eps;
    Poly<T,X> D = this->Deriv();
    T y = this->eval(x);
    int count = 0;
    // eps 1 is temporary
    while ( abs(y) > eps1 ) {
      y = this->eval(x);
      x -= y / D(x);
      // relax convergence
      if ( count++ && !(count % n_relax))
	eps1 *= 2;
      if ( count > n_max ) {
	cerr << "Newton-Raphson convergence problems" << endl;
	cerr << "P(" << x << ") = " << y << ", eps =  " << eps1 << endl;
	break;
      }
    }
    return x;
  }
};

// Legendre Polynomials
template <class T,class X> class Legendre : public Poly<T,X>
{
  using Poly<T,X>::n;
  using Poly<T,X>::C;
public:
  Legendre() : Poly<T,X> () { }
  Legendre(int m) : Poly<T,X>(m) {
    if (!m) {
      C[0] = 1;
    } else if (m == 1) {
      C[0] = 0;
      C[1] = 1;
    } else {
      Legendre<T,X> P2(m-2),P1(m-1); // recursion
      T C2,tmp = 1;
      Poly<T,X> P = P1.Shift(1);
      tmp /= m;
      C2 = tmp - 1;
      P *= 2 - tmp;
      P += P2 * C2;
      for (int k=0; k <= m; ++k)
	C[k] = P.Coeff(k);
    }
  }
  Vector<T> Roots ( T p_eps=0 ) {
    // eps - default convergence criterion
    // convergence criterion will depend on template
    // argument

    T eps = max ( numeric_limits<T>::epsilon(), p_eps );

    // arbitrary precision has an eps of zero so use
    // long double as the lower limit
    double ldeps = numeric_limits<long double>::epsilon();

    if ( eps < ldeps )
      eps = ldeps / 4;

    //cout << "epsilon = " << eps << endl;

    if ( n == 0 ) {
      Vector<T> vr(1);
      return vr;
    } else if ( n == 1 ) {
      Vector<T> vr(1);
      vr[0] = 0;
      return vr;
    } else if ( n == 2 ) {
      Vector<T> vr(2);
      T tmp = 1 / sqrt(3);
      vr[0] = -tmp;
      vr[1] =  tmp;
      return vr;
    } else { // n > 2
      Vector<T> vr(n); // for roots
      T m1=-1;
      // denominator
      Poly<T,X> Den(2);
      Den[2] = 1;
      Den[1] = 0;
      // working poly
      Poly<T,X> P = *this;

      // zero root - odd poly
      if (odd(n))
	vr[n/2] = 0;

      // find roots in pairs
      T x=-1; // start at -1
      for (int j=0; j < n/2; ++j)
	{
	  // find the root
	  if ( P.order() > 2 ) {
	    x = P.NewtonRaphson ( x, 4*eps, 8, 32 );
	  } else if ( P.order() == 2 ) {
	    // quadratic formula
	    x = -P[1] + sqrt(P[1]*P[1] - 4*P[0]*P[2]);
	    x /= P[2] * 2;
	  }

	  //cerr << "root: " << x << endl;
	  // roots are in +/- pairs
	  if (abs(x) < eps) { // zero root - should not get here
	    x = 0;
	  } else {
	    vr[j] = x > 0 ? m1 * x : x;
	    vr[n-j-1] = m1 * vr[j];
	  }
	  if ( P.order() <= 2 )
	    break;

	  Poly<T,X> N = P;
	  // divide poly
	  Den[0] = -x*x;
	  P = N / Den;
	}
      // refine the root calculation
      // roots in pairs
      for (int j=0; j < n/2; ++j) {

	x = NewtonRaphson ( vr[j], eps, 8, 32 );

	vr[j] = x > 0 ? m1 * x : x;
	vr[n-j-1] = m1 * vr[j];
      }      
      return vr;
    }
  }
};

// Tchebycheff Polynomials
template <class T,class X> class Tchebycheff : public Poly<T,X>
{
  using Poly<T,X>::n;
  using Poly<T,X>::C;
public:
  Tchebycheff() : Poly<T,X> () { }
  Tchebycheff(int m) : Poly<T,X>(m) {
    if (!m) {
      C[0] = 1;
    } else if (m == 1) {
      C[0] = 0;
      C[1] = 1;
    } else {
      Tchebycheff<T,X> P2(m-2),P1(m-1); // recursion
      Poly<T,X> P = P1.Shift(1);
      P *= 2;
      P -= P2;
      for (int k=0; k <= m; ++k)
	C[k] = P.Coeff(k);
    }
  }
  Vector<T> Roots ( T p_eps=0 ) {
    // eps - default convergence criterion
    // convergence criterion will depend on template
    // argument

    T eps = max ( numeric_limits<T>::epsilon(), p_eps );

    // arbitrary precision has an eps of zero so use
    // long double as the lower limit
    double ldeps = numeric_limits<long double>::epsilon();

    if ( eps < ldeps )
      eps = ldeps / 4;

    //cout << "epsilon = " << eps << endl;

    if ( n == 0 ) {
      Vector<T> vr(1);
      return vr;
    } else if ( n == 1 ) {
      Vector<T> vr(1);
      vr[0] = 0;
      return vr;
    } else if ( n == 2 ) {
      Vector<T> vr(2);
      T T2=2;
      T tmp = 1 / sqrt(T2);
      vr[0] = -tmp;
      vr[1] =  tmp;
      return vr;
    } else { // n > 2
      Vector<T> vr(n); // for roots
      T m1=-1;
      // denominator
      Poly<T,X> Den(2);
      Den[2] = 1;
      Den[1] = 0;
      // working poly
      Poly<T,X> P = *this;

      // zero root - odd poly
      if (odd(n))
	vr[n/2] = 0;

      // find roots in pairs
      T x=-1; // start at -1
      for (int j=0; j < n/2; ++j)
	{

	  // find the root
	  if ( P.order() > 2 ) {
	    x = P.NewtonRaphson ( x, 4*eps, 8, 32 );
	  } else if ( P.order() == 2 ) {
	    // quadratic formula
	    x = -P[1] + sqrt(P[1]*P[1] - 4*P[0]*P[2]);
	    x /= P[2] * 2;
	  }

	  //cerr << "root: " << x << endl;
	  // roots are in +/- pairs
	  if (abs(x) < eps) { // zero root - should not get here
	    x = 0;
	  } else {
	    vr[j] = x > 0 ? m1 * x : x;
	    vr[n-j-1] = m1 * vr[j];
	  }
	  if ( P.order() <= 2 )
	    break;

	  Poly<T,X> N = P;
	  // divide poly
	  Den[0] = -x*x;
	  P = N / Den;
	}
      // refine the root calculation
      // roots in pairs
      for (int j=0; j < n/2; ++j) {

	x = NewtonRaphson ( vr[j], eps, 8, 32 );

	vr[j] = x > 0 ? m1 * x : x;
	vr[n-j-1] = m1 * vr[j];
      }      
      return vr;
    }
  }
};

// Hermite Polynomials
template <class T,class X> class Hermite : public Poly<T,X>
{
  using Poly<T,X>::n;
  using Poly<T,X>::C;
public:
  Hermite () : Poly<T,X> () { }
  Hermite (int n) : Poly<T,X>(n) {
    if (!n) {
      C[0] = 1;
    } else if (n == 1) {
      C[0] = 0;
      C[1] = 2;
    } else {
      // Hx(x) = 2 x Hn-1(x) - 2 (n-1) Hn-2(x)
      Hermite<T,X> P2(n-2),P1(n-1); // recursion
      Poly<T,X> P = P1.Shift(1);
      P *= 2;
      P2 *= 2 * (n-1); 
      P -= P2;
      for (int k=0; k <= n; ++k)
	C[k] = P.Coeff(k);
    }
  }
  Vector<T> Roots ( T p_eps=0 ) {
    // eps - default convergence criterion
    // convergence criterion will depend on template
    // argument

    T eps = max ( numeric_limits<T>::epsilon(), p_eps );

    // arbitrary precision has an eps of zero so use
    // long double as the lower limit
    double ldeps = numeric_limits<long double>::epsilon();

    if ( eps < ldeps )
      eps = ldeps / 4;

    //cout << "epsilon = " << eps << endl;

    if ( n == 0 ) {
      Vector<T> vr(1);
      return vr;
    } else if ( n == 1 ) {
      Vector<T> vr(1);
      vr[0] = 0;
      return vr;
    } else if ( n == 2 ) {
      Vector<T> vr(2);
      T T2=2;
      T tmp = 1 / sqrt(T2);
      vr[0] = -tmp;
      vr[1] =  tmp;
      return vr;
    } else { // n > 2
      Vector<T> vr(n); // for roots
      T m1=-1;
      // denominator
      Poly<T,X> Den(2);
      Den[2] = 1;
      Den[1] = 0;
      // working poly
      Poly<T,X> P = *this;

      // zero root - odd poly
      if (odd(n))
	vr[n/2] = 0;

      // find roots in pairs
      T x=-10; // start at -2
      for (int j=0; j < n/2; ++j)
	{
	  // find the root
	  if ( P.order() > 2 ) {
	    x = P.NewtonRaphson ( x, 8*eps, 4, 32 );
	  } else if ( P.order() == 2 ) {
	    // quadratic formula
	    x = -P[1] + sqrt(P[1]*P[1] - 4*P[0]*P[2]);
	    x /= P[2] * 2;
	  }

	  //cerr << "root: " << x << endl;
	  // roots are in +/- pairs
	  if (abs(x) < eps) { // zero root - should not get here
	    x = 0;
	  } else {
	    vr[j] = x > 0 ? m1 * x : x;
	    vr[n-j-1] = m1 * vr[j];
	  }
	  if ( P.order() <= 2 )
	    break;

	  Poly<T,X> N = P;
	  // divide poly
	  Den[0] = -x*x;
	  P = N / Den;
	}
      // refine the root calculation
      // roots in pairs
      for (int j=0; j < n/2; ++j) {

	x = NewtonRaphson ( vr[j], eps, 8, 64 );

	vr[j] = x > 0 ? m1 * x : x;
	vr[n-j-1] = m1 * vr[j];
      }      
      return vr;
    }
  }
};

template <class T,class X> ostream &
operator<<( ostream &os, Poly<T,X> &P) {
  P.Write(&os);
  return os;
}

typedef Poly<double,double> DPoly;
typedef Legendre<double,double> LDPoly;
//typedef Poly<double,complex> DCPoly;
typedef Poly<float,float> SPoly;

#endif // _POLY_HPP
